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Abstract 

The temporal component of the heavy-light axial vector current is constructed 
to one-loop order in perturbation theory and to order 1/M, where M is the 
heavy quark mass, in terms of operators suitable for use in lattice simulations 
of B and D mesons. The 0(a)-improved clover action is used for the massless 
light quark, where a is the lattice spacing, and propagation of the heavy quark 
is described by a nonrelativistic lattice action. 
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I. INTRODUCTION 



Studies of B and D meson decays shed light on many important aspects of particle 
physics, such as the extraction of the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements, 
CP violation, rare decays, and hints of physics beyond the standard model. Such decays 
involve not only the electroweak currents, but also strong interaction dynamics in the form of 
hadronic matrix elements. Calculation of these matrix elements is crucial for high precision 
tests of the standard model. 

A very promising approach to determining hadronic matrix elements from first principles 
in QCD is provided by numerical simulations of hadrons using lattice gauge theory. Studies 
of leptonic and semi-leptonic decays using lattice simulations have already been done using 
various techniques to handle heavy quarks on the lattice, such as the static approach [|TJ, 
nonrelativistic QCD (NRQCD) [0,[§, and a reformulation of conventional Wilson and clover 
lattice fermions Q. In each approach, the electroweak currents must be constructed in 
terms of appropriate operators comprised from the gluon and quark fields of the lattice 
theory by calculating the renormalization factors which relate the hadronic matrix elements 
in the lattice- regulated theory to those in the MS scheme. These factors are crucial for 
both decay constant and semi-leptonic form factor calculations. In this paper, we construct 
the time component of the heavy-light axial vector current in terms of operators suitable 
for simulations [|5]|| in which the heavy quark is treated using lattice NRQCD and the light 
quark propagates according to the 0(a)-improved clover action J7|, where a is the lattice 
spacing. The standard Wilson action is used for the gluons. This construction is carried 
out to one-loop order in perturbation theory and to order 1/M, where M is the mass of the 
heavy quark. 

The expansion of the heavy-light axial vector current in terms of appropriate lattice 
operators is achieved by matching relevant scattering amplitudes in perturbation theory. 
Consider the heavy-light pseudoscalar meson decay constant fp$ defined in some continuum 
renormalization scheme, such as MS, by, 

(0\A fl \PS(p)) QCD = ifpsP», (1) 

where is the heavy-light axial vector current and \PS(p)) is a pseudoscalar meson state 
of momentum p. As will be seen, the axial vector current operator can be written as a sum 
of operators J^t m ^ ne lattice-regulated theory. For instance, one finds that the temporal 
component, \x = 0, is a sum of three operators which contribute through 0(a s ,l/M,a), 
where a s is the QCD coupling: 

A = £ C^aM) 4l t + 0(a 2 s ,a 2 ,l/M 2 ,a s a/M). (2) 

3=0,1,2 

Since the role of the dimensionless coefficients Cj is to compensate for low-energy effects in 
the current operator from the loss of short-wavelength QCD modes in the lattice theory, one 
expects that, due to asymptotic freedom, they may be computed to a good approximation in 
perturbation theory. They are fixed by requiring that scattering amplitudes involving these 
current operators agree in continuum and lattice QCD to a given order in a s , a, and 1/M. 
In a perturbative determination, a process involving quarks and gluons as asymptotic states 
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may be used for this matching; one need not consider hadron scatterings which are much 
more complicated to evaluate. Note, however, that the Cj coefficients are independent of the 
process chosen. A one-loop evaluation of the Cj coefficients then involves the following steps: 
(1) select a quark-gluon scattering process induced by the heavy-light current and calculate 
the one-loop amplitude for this process in continuum QCD; (2) expand the amplitude in 
powers of 1/M; (3) identify operators in the lattice theory that reproduce the terms in this 
expansion; (4) calculate the one-loop mixing matrix of these operators in the lattice theory; 
(5) adjust the Cj coefficients to produce a linear combination of lattice current operators 
whose one-loop scattering amplitude agrees with that from the MS current to a given order 
in 1/M and a. Note that the lattice currents are matched directly to the continuum QCD 
currents; we do not need to use a two-step matching procedure in which the lattice currents 
are first matched to those in a continuum heavy-quark effective theory P~CT], which are 
then matched to the full QCD currents [|i~2]1 . 

This paper is organized as follows. First, in Sec. [HJ we discuss various issues relevant 
to our matching procedure, such as the choice of expansion parameter in the perturbation 
series, the definition of the heavy quark mass, the regularization of infrared divergences, and 
the selection of a renormalization scheme. Next, a scattering process appropriate for con- 
structing the heavy-light axial vector current in the lattice NRQCD-clover theory is chosen. 
We then discuss the continuum QCD calculations, and the necessary current operators in the 
lattice regularization scheme are identified. Sec. |T| describes the mixing matrix calculation 
in the lattice theory. The scattering amplitudes in the continuum and lattice schemes are 
matched in Sec. [IV], and our main results, the values of the coefficients Cj for various heavy 
quark mass values, are presented. Renormalon ambiguities are discussed, and subtleties 
with infrared divergences for Wilson light quarks are pointed out. 

II. CONTINUUM CALCULATION AND OPERATOR IDENTIFICATION 

Our goal here is the expansion of the temporal component of the (renormalized) axial 
vector current in terms of (bare) operators suitable for use in lattice simulations of B and D 
mesons. Note that because the axial vector current is partially conserved in QCD, its oper- 
ator does not actually require renormalization. This expansion is determined by matching 
suitable matrix elements evaluated both in continuum QCD and in the lattice theory. Since 
this expansion is an operator relation, the expansion coefficients are independent of the ex- 
ternal states of the matrix elements chosen. Of course, the normalizations of the external 
states used in the continuum and lattice matrix elements must match in order to expose the 
required operator relation. Since the external states in our perturbative calculations are free 
quark states, this can be done at tree level by matching the normalizations of the continuum 
and lattice Dirac spinors (to the appropriate order in a) . The use of on-shell wave function 
renormalization for both the heavy and light quark fields is then the simplest way to ensure 
this matching beyond tree level. 

In matching the 1/M expansions of heavy-light matrix elements in the continuum and 
lattice theories, it is important to establish the relations between the quark mass parameters 
(both heavy and light) in each. Equality of the lattice and continuum mass parameters can 
be ensured in perturbation theory using the pole mass definition. The pole mass is gauge- 
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invariant and, within the context of perturbation theory, is a physical observable. Hence, in 
our calculations, we use the pole mass for M and also for the light quark mass m = 0. In 
using the pole mass definition for our quarks, we compute the continuum and lattice matrix 
elements using on-shell quark mass renormalization. Note that since the pole mass is not 
defined outside of perturbation theory, it is not a useful quantity for lattice simulations; 
the bare mass Mo appearing in the lattice action is much more suitable. Thus, once we 
determine the coefficients Cj(a s ,aM) in terms of M, we re-express our results in terms of 
aM Q using the perturbative relation between M and M . 

We must also know the relationship between the expansion parameters used in the con- 
tinuum and lattice perturbative expressions in order to carry out the operator matching. 
For example, we could use a-jjg in the continuum calculations and a coupling ay defined in 
terms of the momentum-space static quark potential V(q) in the lattice computations. The 
two couplings are related by ay = a jjg + 0(a 2 s ). Our one loop calculations do not make 
use of the 0(a 2 s ) correction in the relation between the lattice and continuum expansion 
parameters; hence, as long as we use couplings which match at leading order, we can simply 
use a s to refer to both the lattice and continuum couplings. Of course, when numerical 
values for the coefficients Cj(a s , aM) are needed, we must specify a renormalization scheme 
for a s . 

The coefficients Cj(a s ,aM) depend on the ultraviolet regulator but contain no infrared 
divergences. However, our choice of on-shell mass and wave function renormalization will 
lead to infrared divergences at intermediate stages in our matrix element calculations. In 
order to demonstrate the cancellation of these divergences in the Cj coefficients, the same 
infrared regulator must be used in the evaluation of both the continuum and lattice matrix 
elements. Since the triple gluon vertex plays no role in any of our Feynman diagrams, we 
introduce a gluon mass A for this purpose. 

Errors in our operator matching are measured in terms of powers of the lattice spacing 
a, the QCD coupling a s , and the inverse heavy quark mass 1/M. However, NRQCD is 
a non-renormalizable effective field theory in which the lattice spacing cannot be taken to 
zero; rather, a = 0(1/M) must be imposed. Hence, our error determinations must take into 
account the fact that aM is of order unity. Since we wish to construct a current operator 
which is correct through 0(1/M), we must remove all 0(a) errors. We also remove all 
0(a/M) errors since it is easy and desirable to do so. The 0(a,l/M,a/M) errors can be 
eliminated at tree-level in perturbation theory. The one-loop calculations are then used to 
remove the 0(a s a) and (a s /M) errors. The remaining 0(a s a/M) errors are much smaller 
than 0(a) and 0(a s a) and are expected to be comparable to 0(a s a 2 ). 

A scattering process which is suitable for our study of the heavy-light axial vector current 
is depicted in Fig. [l|. In this process, an incoming heavy quark | h(p) ) with momentum p 
scatters off a heavy-light current Aq into an outgoing light quark | q(p') ) with momentum 
p' . The axial- vector current operator is given by A^(x) = q(x) 757^ h(x) in terms of the light 
quark field q(x) and the heavy quark field h(x). The first step in our study is the calculation 
of the amplitude for this continuum-QCD process to one-loop order in perturbation theory. 
Using the on-shell mass and wave function renormalization scheme in Feynman gauge and 
expanding in 1 /M (except for the Dirac spinor of the heavy quark) , one finds 
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(q(p')\A \h(p)) 
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where Uh{p) and u q (p') are the standard spinors for the heavy and light quarks, respectively, 
which satisfy the Dirac equation. The light quark mass is set equal to zero. Our conventions 
for the Dirac 7-matrices in Minkowski space are given in the appendix. Ultraviolet diver- 
gences are regulated using dimensional regularization and fully anti-commuting 75 matrices 
are used. As previously mentioned, we use a gluon mass A in order to regulate infrared 
divergences. The three diagrams which contribute to this amplitude at one-loop order are 
shown in Fig. [I]. 

In lattice NRQCD, the heavy quark is described in terms of a two-component (in spin 
space) field ip(x). The Dirac field h(x) is related to ip (and the antiquark field ip) by a 



unitary Foldy-Wouthuysen transformation [13 



h(x) 




(5) 



This transformation decouples the upper and lower components of the Dirac field, thereby 
separating the quark field from the antiquark field. To facilitate the identification of lattice 
NRQCD operators capable of matching Eq. ^, we similarly transform the external state 
spinor Uh{p) (with normalization u\ t Uh = 1) into a nonrelativistic Pauli spinor: 



Uhip) 



2M 



(t-p) 



u Q (p) + 0(1/M 2 ), 



where 




(6) 



(7) 



and Uq is a two-component external state spinor depending only on the spin of the heavy 
quark. Using Eq. |, the relation %uq(p) = uq(p), and the Dirac equation for the light 
quark u q (p')p' = Ug^p')^ ■ p'%, Eq. |^ may be written 
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( q(p') | A) | hip) ) QCD = Vo + V? "I + % A ^2 + 0{al 1/M 2 
where Y = %%, and 

= u q (p') T uq(p), 

n 2 



- m ?(p') m q(p)> 



-^GO 



(8) 

(9) 
(10) 

(11) 



The coefficients in Eq. |8] may be written 
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Having obtained the 1/M expansion of the above scattering amplitude in continuum 
QCD, the next step is to identify operators in the lattice theory which can reproduce the 
terms in this expansion. An inspection of Eq. || suggests immediately that matrix elements 
of the following three lattice operators should be considered: 



^A,lat( x ) 



q{x)TQ(x), 

^(*)r 7 .vQ(*), 
-Lg (a 7 .^rQ(*), 



(14) 
(15) 

(16) 



where q(x) is now the light quark field in the lattice theory, M is the bare heavy quark 
mass, and Q(x) is related to the heavy quark field ip(x) in lattice NRQCD by 



Q(x) 




(17) 



We use the bare quark mass Mo in the above definitions since it is the natural mass to use 
in lattice simulations and because the pole mass M is not well defined beyond perturbation 
theory. The covariant finite difference operators V M and V M are defined as usual in terms of 
the link variables U^x) which are the parallel transport operators from sites x to neighboring 
sites x+a^ in the gauge field. The definitions of these operators are given below, along with 
other derivative operators which will be needed later: 
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(18) 
(19) 

(20) 
(21) 



where 0(x) is an operator defined at lattice site x with appropriate color structure, and Uq 
is the mean link parameter introduced by the tadpole improvement procedure [ 14| . Note 
that these lattice operators are defined in Euclidean space; our Euclidean space conventions 
are outlined in the appendix. 



III. LATTICE CALCULATION 

In this section, we describe the one-loop calculation in the lattice theory of the mixing 
matrix defined by 

( ?(?') I Jfut I %) )iat = E % + 0{al 1/M 2 , a 2 , a s a/M), (22) 



where flj are the Euclidean-space counterparts of the flj defined in Eqs. |9|-|11|. First, the 
lattice actions used in these calculations are specified. The necessary Feynman diagrams are 
then presented, and their evaluation in lattice perturbation theory is outlined. 
For the heavy quark, we use the following NRQCD action density 0: 

aC-NRQCD = ^(x) lf)(x) 

t , J a5H\ ( aH \ n uUx) ( aH \ n ( a5H\ , . 



where 



A( 2 ) 

F » = -2A4' (24) 

W =- C "4 CT ' R (28) 

The positive integer n is introduced to stabilize the highest momentum modes in the heavy 
quark propagator 0; the condition, n > 3/aM , has proven to be a reliable guide. Note 
that the above NRQCD action does not include the heavy quark mass term. The QCD 
coupling g is related to a s in the usual manner, a s = g 2 /(47r), and <x,- are the standard Pauli 
spin matrices. At tree level, cb = 1; the one-loop contribution to cb is an 0(a 2 ) effect in 
our mixing matrix calculation and hence can be ignored here. The chromomagnetic field 
is given by Bj(x) = — \eji m Fi m (x), where the Hermitian and traceless field strength tensor 
Fpvix) is defined at the sites of the lattice in terms of clover-leaf operators: 
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= T~4 H C4(x)L^(x+a a )t/_ a (x+a Q! +a (3 )C/_^(a:+a ( g), (26) 

with {(a, (3)}^ = {(jjl, u), (is, -//), (-//, -1/), (-1/, //)} for // ^ 1/. 
For the light quarks, we use the clover action 0, 

aCught = qfq - a^-qV^q + m qq - iga^Yl ( 27 ) 

where ^ = X)„ 7 M V^, m is the bare light quark mass, a^ v = \ [7^, 7„], and we set the Wilson 
parameter r — 1. The one-loop correction to the clover coefficient is an 0(a 2 s ) effect in our 
matching calculation and can be neglected. 

Lattice perturbation theory calculations are much more laborious than those in contin- 
uum perturbation theory. Not only do more diagrams contribute to a given process, but 
the complicated functional forms of the propagators and vertex functions (presented in the 
appendix) necessitate the use of numerical methods in evaluating the integrals over internal 
loop momenta. For the scattering process considered here, there are five additional one-loop 
Feynman diagrams in the lattice theory: three are vertex corrections, shown in Fig. |2|, and 
two are extra external leg corrections, shown in Fig. [| The Feynman rules are determined 
by expanding the total lattice action in terms of g using 



UJx) = exp 



iagG,, {x + Y 



(28) 



and u = 1 — a s u^ + 0(al), then Fourier transforming into momentum space. G^(x) is the 
lattice gluon field defined at the midpoints of the links connecting neighboring sites. 

We calculate the amplitudes corresponding to the Feynman diagrams of Figs. [l]-|3| in two 
very different ways and verify that the results agree. In the first method, all spin matrix 
manipulations and derivatives with respect to external momenta are done by hand. The 
resulting integrals are also simplified by hand, including the subtraction of terms to remove 
infrared divergences. The remaining infrared-finite integrals are then done numerically using 
Monte Carlo and adaptive Gaussian quadrature techniques. 

In the second method, nearly all aspects of the calculation are automated. First, the 
propagators and vertex functions are expressed as functions in C++ whose arguments are the 
appropriate four momenta implemented using a class f ourvec. Spin algebra is accomplished 
using explicit matrix representations; this is done by defining new class structures in C++, 
such as pauli and dirac, and overloading all of the necessary arithmetic operators for ease of 
use. Derivatives with respect to external momenta are taken using automatic differentiation 



I5fl . This is implemented by defining in C++ a class Tcomplex which carries out the 
multivariate Taylor series expansions. A variable of type Tcomplex is treated by the end 
user just as if it were a regular complex scalar variable since all appropriate arithmetic 
operators and mathematical functions are overloaded. However, a Tcomplex variable is 
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actually an array containing the function value and all of its derivatives up to some order; 
an indexing member function is used to return a specific term in the Taylor series expansion. 
All derivatives are taken automatically using analytical techniques; the end user never needs 
to write subroutines for taking these derivatives. 

The above tools allow one to easily write C++ functions to compute the integrands for 
all of our one-loop Feynman diagrams. The integrals corresponding to the graphs in Figs. |2| 
and [5] are easily evaluated using Monte Carlo techniques. We use the integration routine 
VEGAS Note that, in lattice perturbation theory, each internal loop momentum k is 

restricted to the first Brillouin zone — 7r < ak^ < n. However, due to infrared divergences, 
the other integrals cannot be directly calculated. We evaluate these integrals in a sequence 
of steps. First, we subtract from the lattice integrand Ii a t{k), where k is the loop momentum, 
the analogous integrand from the continuum theory, I con (k). The integral of this difference 
over the first Brillouin zone is infrared finite and can be evaluated using the Monte Carlo 
method. We then multiply the continuum integrand by a factor f(k) to render the integral 
insensitive to the shape of the first Brillouin zone. This factor must also be chosen such 
that the integral of [1— f(k)] I con (fc) is infrared finite. For example, f(k) = exp(— pk 2 ) where 
p = 7r/2 is often a good choice. We then use the Monte Carlo method to evaluate the 
integral of [1 — f(k)] I con (k). The infrared divergence has then been isolated in the integral 
of f(k) I con (k). Since f(k) is chosen to remove any appreciable sensitivity of the integral to 
the edge of the Brillouin zone, we can change to hyperspherical coordinates and integrate 
over the interior of an infinitely large sphere to a good approximation. This simplifies the 
calculation and allows us to easily identify and analytically manipulate the infrared divergent 
pieces of the integrand. 

Using these two methods, we obtain results for the mixing matrix in Eq. [22], again 
using the on-shell renormalization scheme in Feynman gauge. Note that the nonrelativistic 
external-state spinors ug(p) are identical in the lattice and continuum theories. For the 
clover action, the external light-quark spinors u q (p') differ only at 0(a 2 ). At one-loop, the 
mixing matrix elements may be written 



Zij = da + a s 



1 „ 

-(C<? + Cq) $ij + Cm ^il^j'l + Cij 



(29) 



where C q and Cq are the contributions from the light- and heavy-quark external leg correc- 
tions (that is, from wave function renormalization factors), and Qj denote the contributions 
from the vertex corrections. Our use of an on-shell renormalization scheme with lattice 
operators defined in terms of the bare mass M is responsible for the term proportional to 
C m , where 

M =[l + a s C m ]M + O(a 2 ). (30) 

(i) 

Note that although the current operators J A Ut are defined using M , the pole mass M must 
appear in Qj. The factors in Eq. [29] may be further decomposed: 

4 

Co = Co InaA, 

Q Q 3vr 
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— i ^ m i 

Gi-Cii-CS 1 -^: (31) 

where C q , Cq, C m , and Qj are infrared finite and independent of the tadpole improvement 
factor m , and (jf- and contain the infrared divergences and tadpole improvement con- 
tributions, respectively, from the vertex corrections. Contributions to C q and C m from the 
tadpole improvement counterterms are denoted by Cj l and C^ 1 , respectively. 



IV. RESULTS 

To complete the operator matching, we transform Eq. M from Minkowski to Euclidean 



space, use Eqs. £2] and g9], then peel off the dependence on the external states. This yields 
the following operator relation: 
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+ 0(a,,a^,l/M ,a s a/M). 



(32) 



As expected, we find that the infrared divergences from the various terms in the expansion 
coefficients cancel. We now have the desired expansion coefficients Cj of Eq. ^| in terms of 
the pole mass M. As previously mentioned, the bare mass Mo is a much more convenient 
mass parameter. Our results are easily expressed in terms of aMo using Eq. which, in 
this case, simply amounts to replacing M by M : 
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1 + a s 
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(C q + Cf + Cq) - T 



J 



(0) 
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+ 0(a s , a , 1/M 2 ,a s a/M), 



(33) 



where B = ln(aM )/7r - 1/(4tt), B 1 = ln(aM )/7r - 19/(12tt), B 2 = 4/vr, r = Coo + C10, 
T i — Coi + Cii) an d t 2 = Co2 + Ci2- The factors arising from tadpole improvement counterterms 
are given by 



C 



c 



TI 
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TI 

m 

TI 



-U 



(2) 
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-v {2) 



2naM 
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(2) 



(34) 



For the usual plaquette definition u = (^Trf/^) 1 / 4 in the Wilson gluonic action, u = tt/3. 
For massless clover quarks, C q = 1.030. Results for r , Ti, r 2 , Cq and C m for various values 
of aMo are listed in Table |. 
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Explicit calculation reveals that r 2 behaves as —2aMo(disc as Mo becomes large, where 
(disc is found numerically to be (disc — 1-00(1). This factor can be viewed as arising from a 
discretization correction J^ut to ^ii at , 

= aq(x) 1 -^TQ(x), (35) 
analogous to the 0(a) correction a d^P, where P is the pseudoscalar density, to the axial 



current in light-light quark systems JT7j. We can then define an improved current operator 

7 (0)imp _ j(0) n j{disc) /o fi \ 

where we write 



c -=n l+ ^J' (37) 

This decomposition is not unique, and hence, we leave (a as a free parameter. Different 
choices of (a lead to different 0(a 2 ) contributions being included in Eq. [33]. Taking this into 
consideration, Eq. |33] can be written 

A = (1 + a s p ) J { A °}l™ p + (1 + a 8 fn) Jam + ^{p2 - (a) Jam 

+ 0(a 2 s ,a 2 ,l/M 2 ,a s a/M), (38) 

so that the Cj coefficients of Eq. become Co = l+a s po, C\ = l+a s pi, and C2 = a s (p2 — (a)- 
This equation is our final result. Note that at 0(a s ), the (a dependence in Ja^ P cancels 

that in C 2 JaLk,- Numerical results for p Q , p 1; and p 2 are given in Table [TT]. 

In Ref. 0, B meson decay constants were computed in simulations in which the heavy 
quark propagation was described by an NRQCD action different from Eq. in particular, 
higher order corrections were included. The values of po, pi, and p2 for the action used in 



these simulations are given in Table [III . 



In order to use Eq. |38] in a simulation, a value for a s must be specified. To do this, 
one must first choose a renormalization scheme for a s ; in the case of a running coupling, 
a means of setting the scale must then be devised; lastly, the value of the coupling at 
some reference scale must be determined. A coupling a>y(q*) defined in terms of the short- 
distance static potential with a scale-setting prescription based on the mean value theorem 
is advocated in Ref. jnj]; the value of ay(3.402/a) can be obtained from measurements 



of the average plaquette. Unfortunately, we have not computed the q* scales for the p , 



Pi, and P2 coefficients. However, based on findings in Refs. and |Tj|, we expect that 
q* ~ 2/ a. An alternative choice for the expansion parameter is a non-running boosted 
coupling ab = 6/(47T/9mo); the value of this coupling is typically comparable to oty(7r/a). 

For the range of aM values considered here, our results show no evidence that pertur- 
bation theory is failing. From Tables |TT] and [TTT], one sees that the one-loop corrections to 
Cq are very small. Values for po are typically near —0.3 for all aMo considered; multiplying 
by a s ~ 0.2 yields one-loop corrections to Co which are approximately 5% of the tree-level 
contribution. For the large range of aM values studied, p\ varies between —0.7 and 0.5; 
hence, one-loop corrections to C\ are never larger than about 15%. The magnitudes of the 
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values for p 2 are significantly larger than those for p and p\. If we set (a = 0, then Ca = oc s 
and C*2 = <y s p2 which, using a s ~ 0.2, varies between 0.2 and 0.9 for the range of aM values 
studied here. Since Ci has no tree-level contribution, at least a two-loop result is needed to 
check the convergence of its perturbative expansion. 

As in all applications of perturbative QCD, nonperturbative contributions to Cj are pos- 
sible. For example, the perturbative expansions of the Cj coefficients contain renormalon 
ambiguities arising from a slight mismatch between the infrared physics of the lattice theory 
and that of continuum QCD. This mismatch is caused by discretization artifacts and the 
truncation of the 1/M expansion. In our calculation, contributions coming from momenta 
small relative to a given infrared scale q <C l/a,M are all suppressed by a 2 g 2 or q^/M 2 
since the infrared structure of our lattice theory is designed to agree with continuum QCD 
through 0(a,l/M). This infrared suppression is not true diagram-by-diagram (some dia- 
grams are even infrared divergent), but the infrared contribution is suppressed by a 2 or 1/M 2 
when all diagrams are summed. For example, in some analyses [HJ, coefficients like Cq are 



separated into the contribution c from the leading operator (here «/i°j at ) and contributions 

Co from mixings with higher-dimension operators (here J^Lt anc ^ ^Aiat)- m our analysis, 
such a separation is artificial and would be a mistake since it induces 0(a, 1/M) ambiguities 
in the separate pieces c and c . These ambiguities cancel when the pieces are recombined, 
leaving contributions to Cq suppressed by a 2 g 2 or q\/M 2 . Thus, any renormalon ambiguities 
are suppressed by the same factors and are at worst comparable to the other truncation 
errors in the analysis. 

Our calculation differs from conventional lattice calculations in that the a 2 truncation 
errors, perturbative and nonperturbative, cannot be made arbitrarily small by reducing 
the lattice spacing [fjj. This is because NRQCD has nonrenormalizable interactions whose 
couplings do not vanish as a — > 0. The perturbative expansions for Cj have power-law terms 
of the form a s /(aM) n which ruin the convergence of perturbation theory if a is taken to 
zero. In practice, this problem is avoided by ensuring that l/(oM) is not large. Our inability 
to take a —>■ in NRQCD fundamentally limits the precision of our 0(a,l/M) accurate 
formalism. If improved precision is needed, we must reduce the truncation errors by using 
more accurate discretizations of the lattice currents and action. In this way truncation 
errors, renormalon ambiguities, etc. are pushed off to 0(a 3 , 1/M 2 ) or higher. 

This approach differs dramatically from that used in the static-quark expansion; in the 
latter approach, problems with power divergences arise when the limit a — > is attempted 
in a brute-force manner. The role of the power-law terms in the Cj coefficients is to cancel 
unphysical power-divergent contributions in the lattice NRQCD current matrix elements, 
thereby ensuring the correct matching to full continuum QCD through the relevant order 
in a s and 1/M. If the unphysical terms behaving as a s /(aM) n in the Cj coefficients and 
the current matrix elements become comparable to or larger than the physical terms of 
0(Aqcd/M), this cancellation becomes a very delicate issue and the need for nonperturba- 
tive subtraction methods may arise. Such is the case in the static-quark expansion; taking 
the brute-force limit a — > obscures the physically relevant contributions by amplifying the 
power-divergent unphysical terms. Here, this problem is circumvented by keeping the lattice 
spacing large enough that l/(oM) never becomes large. The success of low-order pertur- 
bation theory in cancelling the unphysical terms depends, of course, on the relative size of 
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the physical contributions. Recent numerical simulations find no evidence of large 

power-law contamination in the current matrix elements for the range of aM appropriate 
for 6-quark physics, suggesting that power-divergent nonperturbative and higher-order per- 
turbative contributions to the matching coefficients are also not large. This finding applies 
to the temporal component of the heavy-light axial vector current. For other operators, the 
situation could be different and each case must be examined separately. 

The results in Table [TT] are suitable for values of aMo appropriate for current simulations 
of heavy-light systems using NRQCD. However, terms proportional to In (aMo) cause the 
0(a s ) contributions to C and C\ to become large as aM becomes large. In such cases, 
the renormalization group should be used to improve upon the estimates from one-loop 
perturbation theory. Since the left-hand side of Eq. || is independent of the lattice spacing, 
it follows that, neglecting 0(a 2 ) terms, 



d_ 

da 



Y,C k (a s ,aM ) (jjg at > 



L k 



0. 



(39) 



Using Eq. ^ and dQj/da = 0, a renormalization group equation for the Cj coefficients, 
collected into a vector C, can be obtained which describes the change in Cj as the lattice 
spacing is varied: 



'a± + 7 tr ) C = 0, 
. da 



where the anomalous dimension matrix is given by 

d 7 " 

da 



7ij(a s ,aM ) 



7- 1 



(40) 



(41) 



In the limit of large aMo, we find 7 = diag(— a s /ir, —a s /ir, —a s /ir) and Eq. f|0] can be 
easily solved. First, express the Cj coefficients as a function of a s (a) and aM instead of a s (a) 
and aMo (a) since M is fixed and does not run with a. If the lattice theory is renormalized in 
such a way that the renormalization group /3-function /3rg — ~ a da s {a^ 1 ) / da is independent 
of M and mo, then 



PrgM = -2/3 a 2 s - 2/3x a s 3 + 0{a% 



(42) 



where (3 = (11 - §iV/)/(47r) and ^ = (102 - f A^ / )/(16vr 2 ), for N f light quark flavors. 
For large, fixed M, Eq. [H] then tells us that the change in Cq and C\ in going from an 
initial lattice spacing a\ to a smaller lattice spacing a 2 is given, in leading logarithmic 
approximation, by 



q(q 2 M) 
Ci{ ai M) 



a B (l/o2) 
a s (l/ai) 



l/(2/3 vr) 



(43) 



for % 



0, 1, where a s (fi) is the familiar QCD running coupling; its two- loop form is 

2 / A 2 \ , Pi 



a 3 (ji) 



PoMv 2 /^ 



+ -ilnln( y u 2 /A 2 ) 

Po 



(44) 
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Let a 2 M = oMq and a\M = aMo, then 

a s (l/a 2 ) 



a s (l/ai) 



1 - a s /5 In 



« S (M 



7T + C(a 



a s (M ) 

Hence, to the order at which we are working, Eq. |43| is equivalent to 



(45) 



Cj{aM' Q ) 
a (aM ) 



a s (M ) 



-l/(2/3„7r) 



(46) 



for large oM and aM' Q . When aM is not large, the dependence of 7^ on aM Q cannot be 
neglected and Eq. |46| must be suitably modified. 

As an aside, we mention the following fact. If one wishes to use the Wilson action instead 
of the clover action for the light quark, one cannot include the discretization correction term 



^Ailt ■ If one includes J^iat'' wnne using Wilson light quarks, one finds an uncancelled loga- 
rithmic infrared divergence upon attempting to match the continuum and lattice scattering 
amplitudes. This divergence is removed by including contributions from the 0(a)-correction 
term in the clover action. 



Adisc) 



V. SUMMARY 

In this paper, the temporal component of the heavy-light axial vector current was 
expanded in terms of lattice operators suitable for use in simulations of B and D mesons. 
The expansion was carried out to 0(1/M) by matching relevant scattering amplitudes to 
one-loop order in perturbation theory. The (massless) light quark was described in the 
lattice theory using the 0(a)-improved clover action of Eq. |27], and the NRQCD action of 
Eq. |23] was used to treat the heavy quark. The standard Wilson action was used for the 
lattice gluons. The expansion of the heavy-light current was found to be 

A, = C Ji°jr + Ck J&t + C2 jgL + Otf, a 2 , 1/M 2 , a s a/M), (47) 

where the lattice operators are defined in Eqs. [14], [16], and |36|, and C = 1 + a s Po, 
Ci = 1 + a s pi, C 2 = a s (p 2 - CO and C A = ot s [l + ( A /(2aM )]. Values for p , Pi, and p 2 
are listed in Table [□] for various bare heavy quark masses, (a remains as a free parameter; 
different choices for the value of (a lead to different 0(a 2 s ) contributions being included 
in Eq. |17]. The one-loop corrections for Cq and C\ were shown to be small relative to 
the tree-level contributions for the range of aM studied; no evidence of a breakdown in 
perturbation theory was found. Since C 2 has no tree-level contribution, at least a two-loop 
calculation would be needed to check the behavior of its perturbative expansion. Renormalon 
ambiguities were argued to be at worst of the same order as our other systematic errors. 
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Our results have already been used in B- meson simulations || measuring More 
recently, B meson decay constants were computed || in simulations using an NRQCD 
action with higher-order interactions not included in Eq. the values of p , pi, and p 2 for 
the action used in these simulations are given in Table |TJ. In the future, we plan to apply 
the methods described in this paper to the expansion of other currents, such as the vector 
current, in terms of appropriate lattice operators. 
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APPENDIX A: FEYNMAN RULES 



The NRQCD and light quark actions are given in Eqs. ^ and [27] in Sec. |T|. The 
Feynman rules of perturbation theory can be derived from these actions using, for instance, 
the methods of Ref. [[H]]. Most of the rules relevant for the calculations of this article are 




collected in this appendix. Various conventions used in our computations are also outlined. 
To simplify notation, we set the lattice spacing a = 1 in this appendix. 

Minkowski-space quantities are indicated either by a caret, such as 7, or by a subscript 
or superscript (M), such as Xj M \ The metric tensor in Minkowski space is taken to be 
g^ u = diag(l, — 1, — 1, — 1) and the Dirac matrices satisfy {7^,7^} = 2g IJjU . We use the 
Dirac-Pauli representation: 

7° = 7o = 
f = _ 7j = 

where <x,- are the standard Pauli spin matrices. Also, we define 75 = ry ^ 1 ^ 2 ^ 3 . 

Euclidean-space four- vectors are defined in terms of Minkowski-space four-vectors using 
xq = x° = ix Q , M ^ and Xj = x- 7 = xl M s = —Xj , for j = 1,2,3. For the derivative operator, 

d = d° = —id^ and dj = d j = dj = —d^ M y Note that the gauge field G^ M ^ = 

(<P(m),G(m)) Wick-rotates into Euclidean space as a covariant vector, as does the gauge- 
covariant derivative D? M y The Euclidean-space Dirac matrices satisfy {7^,7^} = 2 5^ u and 
are related to their Minkowski-space counterparts by 7o = 7° = 7°, Jj = 7 J = —i^ = ijj, 
and 75 = 7o7i7273 = %■ Also, a UfM = | [7^, 7J. To be consistent with our conventions 
for Euclidean-space 7 matrices, we define Euclidean-space quark-bilinear axial and vector 
currents by A = Aq , = —iA*i M \, Vq = V$ M \ and V j = —iV? M y The chromoelectric 
and chromomagnetic fields are defined in terms of the field strength tensor by Ej = Foj = 
-iF { f } = %E) M) and Bj = -~e lmj Fi m = -~s lmj F^ = B^ M \ where e ijk is the fully 
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antisymmetric Levi-Civita tensor and the field strength tensor is given on the lattice by 

The heavy quark propagator is diagonal in both spin and color, and is given in momentum 
space by 

G Q (k) = (l-e- lk0 F 2n (k)y\ 

The light quark propagator is diagonal in color: 

G q {k) = ii ^2 In sin K + 2r E sin 2 (|/^) + 



m 



In this paper, we work only with the m = and r = 1 case. The gluon propagator is 
diagonal both in color and the Lorentz indices, and is given in Feynman gauge by 

G G (k) = (4^sin 2 (iA; M ) + A 2 ) . 

We now list the vertex factors associated with the interaction of a single gluon with a 
heavy-quark line. Let k! be the outgoing heavy-quark momentum, k be the incoming quark 
momentum, and fi be the polarization index of the emitted gluon. These vertex factor all 
have a color factor of T£ c , where a is the color index of the gluon, and b and c are the color 
indices of the outgoing and incoming quarks, respectively. For ji = 0, the vertex factor is 

-^e-* (fc ' +fc)o/2 F n (k) F n (k'). 

For /i = j = 1, 2, 3, the vertex factor from the p 2 /2M term is 

1 



sin[i(A; , + k)j] [e- ik '«F n (k') + e- ik °F n (k)\ S n {k', k), 
S n (k',k) = E?=oF l (k')F n - l -\k). 



9 2nM u '" l2r " '■' 



For = j from the a ■ B term, the vertex factor is given by 

9 E (75707.) e rsj sm(k' - k) s cos[\{k' - k)j] [e~<F 2n (k') + e~^F 2n {k) 

r,s 

The single-gluon vertex factor for the light quark (clover action) is 
-g cos[§(A/ + k)p] +rsin[\(k' + k)^] + ^E "^ sin ( k ' ~ k) v cos{\{k' - fc) M ]|, 

for all gluon polarizations /i, where k' (k) is the outgoing (incoming) momentum of the light 
quark. 
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Next, the vertex factors associated with the interaction of two gluons with a heavy-quark 
line are given. Let k' be the outgoing heavy-quark momentum, k be the incoming quark 
momentum, and qi, q 2 be the outgoing momenta of the emitted gluons having polarization 
indices fj,±, fi 2 , respectively. Only the factors for jii = jj, 2 (and q± + q 2 = 0) are given, which 
is all that is required for the one-loop tadpole graphs. The color factor for these cases is 
{T a T a ) bc , where a, a' are the color indices of the gluons, and b and c are the color indices 
of the outgoing and incoming quarks, respectively. For jii = fj, 2 = 0, the vertex factor is 

a 2 

_y_ e -i(fc'+fc) /2 pn^ pn^y 

2 

For iii — /j.2 — j, the vertex factor from the p 2 /2M term in the NRQCD action is 

(- 1 



g 



e - ik 'oF n (k') + e- iko F n (k) 



cosg(fc' + fcy S n (k', k) 



+ 



(2nMy 



\4nM 

sin[\(k'+k-q 2 )j] sin[l(k' + k + q 1 ) j ] S g , n (k',k' + q±,k) 



+ 9' 



where 



(2nMf 



sm[l(k' + k-q 2 ) j } sm^k'+k+q^j] S n (k', k - q 2 ) S n {k' + q u k) e ~^ k ' + ^\ 



n-2 



S g ,n(k', k' + qi ,k) = J2 F\k') S rM (k' + Ql , k). 

l>0 

The vertex factor for fi\ = fi 2 = j from the cr • B term is 

n 1 



9 



(4My 



Er,s,r>,s> e-W + ^° F 2n (k' + qi ) ( 75 7 07r ) (75707,0 ersj e r 



x sin(/c' — k + q 2 ) s sm(k' — k + qi) s > cos[^(k' — k + q 2 )j] cos[|(fc' — k + q^j] 
+ terms that vanish for qi + q 2 = 0. 

(i) 

The heavy-light currents, J A lat , are not part of the action and we list their Feynman rules 
separately. Again, let k! and k be the outgoing and incoming quark momenta, respectively, 
and let q^ and /x, denote the momenta and polarizations, respectively, of the emitted gluons. 
Also, let q ext be the momentum carried off by the external heavy-light current, where k = 
k' + qi + qext- We use the notation [O]^ to indicate the vertex factor for n-gluon emission 
from operator O. At tree level, 



J 



(i) 



A, lat 



J 



(2) 



A, lat 



(0) 
(0) 



2M 



2M 



Y^sinWik' + k + q ext ) j \ (r 7j ), 
j 

Y,MW + k- qext )j\ ( 7j r). 



For one gluon emission (/i = j, color factor = T b a c ): 
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J 



J 



(1) ■ 

A,lat 

(2) 
A,lat 



(1) 
(1) 



~ 9 0»7 COS k( k ' + k + Qext)j] (T7j) 



= - 9 



2M 

i 

2M 



cos[ 1 (A; / + k — q e xt)j] (lj r) . 



For two gluon emission (fii — ji 2 — j, color factor = (T a T a ') bc ): 



J 



(i) 



J 



A,lat 
(2) 



(2) 



A.lat 



9 AM ~~ 12 

(2) , i 



sin^A;' + fc + q ex t)j] (Flj) , 



= g2 AM Sin ^^' + fc ~ 9ert ^ irij r) ' 
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FIGURES 




FIG. 1. Feynman diagrams in continuum and lattice perturbation theory which contribute 
at one-loop order to the annihilation by the axial current (cross inside a circle) of an incoming 
heavy quark (double solid line) and the creation of an outgoing light quark (single solid line). The 
exchange of a gluon is denoted by a curly line. 
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FIG. 2. Additional vertex correction diagrams which contribute in lattice perturbation theory 
to the same process as in Fig. [|. 



21 




22 



TABLES 



TABLE I. Values of the coefficients To, t%, T2, Cq, and C m appearing in Eq. 33 for various 
values of the bare heavy-quark mass 0M0 and NRQCD stability parameter n. Uncertainties in the 
determinations of these parameters due to the use of Monte Carlo integration are included. 



aM 


n t c 




T~2 


Cq 




10.0 


1 (J.o2<32(l) — 1.242(4) 


-14.74(2) 


0.2719(8) 


0.872(8) 


7.0 


1 0.7929(1) — 1.266(1) 


-8.82(2) 


0.1847(6) 


0.862(6) 


4.0 


1 0.7290(1) -1.283(2) 


-3.202(8) 


-0.0287(6) 


0.883(6) 


4.0 


2 0.7397(1) -1.298(2) 


-3.476(8) 


-0.0030(6) 


1.087(6) 


3.5 


2 0.7239(1) -1.307(2) 


-2.618(8) 


-0.0651(6) 


1.120(6) 


3.0 


2 0.7052(1) -1.315(2) 


-1.790(8) 


-0.1463(6) 


1.142(6) 


2.7 


2 0.6923(1) -1.324(2) 


-1.314(4) 


-0.2077(6) 


1.172(6) 


2.5 


2 0.6829(1) -1.332(2) 


-1.010(4) 


-0.2562(6) 


1.195(6) 


2.0 


2 0.6576(1) -1.352(1) 


-0.304(3) 


-0.4156(6) 


1.235(4) 


1.7 


2 0.6422(1) -1.371(2) 


0.081(3) 


-0.5503(8) 


1.277(4) 


i a 
1.0 


2 0.6377(1) -1.377(1) 


0.202(3) 


-0.6038(8) 


l.l\)6yo ) 


l.Z 


3 0.6307(1) -1.424(1) 


0.538(2) 


-0.868(1) 


1.00 ( ^DJ 


l.U 


4 0.6365(1) -1.457(1) 


0.683(2) 


-1.074(1) 


1. ( oo{u ) 


n s 
U.o 


5 0.6617(1) -1.502(1) 


0.818(1) 


-1.394(1) 


i.yo4^D j 




TABLE II. Values for the coefficients po, 


pi, and p2 appearing in Eq. 38. 




aM 


n 


Po 


Pi 




P2 


10.0 


1 


-0.2972(4) 


0.315(9) 




-3.982(24) 


7.0 


1 


-0.3368(3) 


0.212(6) 




-3.902(16) 


4.0 


1 


-0.3443(3) 


-0.032(6) 




-3.525(8) 


4.0 


2 


-0.3679(3) 


-0.038(6) 




-3.251(8) 


3.5 


2 


-0.3635(3) 


-0.102(6) 




-3.109(8) 


3.0 


2 


-0.3533(3) 


-0.161(6) 




-2.937(8) 


2.7 


2 


-0.3433(3) 


-0.214(6) 




-2.813(4) 


2.5 


2 


-0.3341(3) 


-0.253(6) 




-2.717(4) 


2.0 


2 


-0.3002(3) 


-0.343(4) 




-2.423(3) 


1.7 


2 


-0.2691(4) 


-0.420(4) 




-2.208(3) 


1.6 


2 


-0.2571(4) 


-0.451(6) 




-2.128(3) 


1.2 


3 


-0.2096(6) 


-0.572(6) 




-1.664(2) 


1.0 


4 


-0.1703(7) 


-0.627(6) 




-1.410(2) 


0.8 


5 


-0.1069(7) 


-0.694(6) 




-1.144(1) 
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TABLE III. Values for the coefficients po, p%, and p2 using the NRQCD action of Ref. Q. 



aM 


n 


Po 


Pi 


P2 


20.0 


1 


-0.1532(5) 


0.69(2) 


-4.85(4) 


12.5 


1 


-0.2432(4) 


0.51(1) 


-4.82(4) 


inn 
1U.U 


i 
i 


n 9779 f /A 
— U.Z 1 1 ^{Q) 


U. 40 1^1 ) 


A 7lV'} N l 
— 4. 1 U[o ) 


7.0 


1 


-0.3174(4) 


0.322(7) 


-4.38(2) 


5.0 


2 


-0.3405(4) 


0.208(7) 


-3.89(1) 


4.0 


1 


-0.3372(3) 


0.140(7) 


-3.643(8) 


2.7 


2 


-0.3375(4) 


0.027(4) 


-2.859(4) 


2.0 


2 


-0.3145(4) 


-0.037(4) 


-2.339(4) 


1.6 


2 


-0.2844(4) 


-0.054(4) 


-1.986(3) 
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